##################################################################################################################################################################
##################################################################################################################################################################
#   PNAS
#   The Effect of Big-City News in Rural America During the COVID-19 Pandemic 
#   Eunji Kim, Michael Shephard, Josh Clinton, Vanderbilt University
#   August 2020
##################################################################################################################################################################
##################################################################################################################################################################
rm(list = ls())
directory <- getwd()

library(stargazer)
library(cobalt)
library(backports)
library(stargazer)
library(ggplot2)
library(ggridges)
library(scales)
library(gridExtra)
library(egg)
library(ggpubr)
library(grid)
library(reshape2)


load(file="full2_PNAS.Rdata")
load(file="full_FINAL.Rdata")
load(file="lucidfull.Rdata")
ind <- lucid 
# Subset to Local News Consumers
datlocal <- subset(ind,ind$MORELOCALNEWS > 1)
# Subset to Non local only
datnolocal <- subset(ind,ind$MEDIA2 == 1 | ind$MEDIA2 == 4)



##############################################################################
#        Main Figure 1 
#############################################################################




f1 <- ggplot(full2[which(full2$top25factor!="NA"),], 
             aes(x = log_dmacases, 
                 y = top25factor , 
                 fill = top25factor),  scale_fill_manual(values =c("#1D3752", "#50BFC3"))) +  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.3,
    aes(fill= top25factor)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) + 
  xlab("Log(COVID-19 cases in DMA County Core)") + ylab("")    +  
  scale_x_continuous(breaks = c(3, 6, 9, 12)) +
  scale_fill_cyclical(values = c("#214D72", "#50BFC3")) +  coord_cartesian(xlim =c(0, 12)) 



f2 <- ggplot(full2[which(full2$top25factor!="NA"),], 
             aes(x = log_casesnew, 
                 y = top25factor , 
                 fill = top25factor),  scale_fill_manual(values =c("#1D3752", "#50BFC3"))) +  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.7,
    aes(fill= top25factor)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) + 
  xlab("Log(COVID-19 cases)") + ylab("")    +  
  scale_x_continuous(breaks = c(3, 6, 9, 12)) +
  scale_fill_cyclical(values = c("#214D72", "#50BFC3")) +  coord_cartesian(xlim =c(0, 12)) 



jpeg("Figures/Figure1.jpeg", units="in", width=14, height=7, res=300)

pdf("Figures/Figure1.pdf")

egg::ggarrange(f2, f1, ncol = 1, heights = c(0.89, 1))

dev.off()


# per capita graph 


f3 <- ggplot(full2[which(full2$top25factor!="NA"),], 
             aes(x = log_dmacasespercapita, 
                 y = top25factor , 
                 fill = top25factor),  scale_fill_manual(values =c("#1D3752", "#50BFC3"))) +  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.3,
    aes(fill= top25factor)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) + 
  xlab("Log(COVID-19 cases in DMA County Core per capita)") + ylab("")    +  
  scale_fill_cyclical(values = c("#214D72", "#50BFC3")) + 
  scale_x_continuous(breaks = c(-12, -9, -6, -3, 0))  + coord_cartesian(xlim =c(-12, 0)) 

full2$log_covidpercapitanew <- full2$log_covidpercapita
full2$log_covidpercapitanew[full2$log_covidpercapita==0] <-  min(full2$log_covidpercapita, na.rm=T)

f4 <- ggplot(full2[which(full2$top25factor!="NA"),], 
             aes(x = log_covidpercapitanew, 
                 y = top25factor , 
                 fill = top25factor),  scale_fill_manual(values =c("#1D3752", "#50BFC3"))) +  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.7,
    aes(fill= top25factor)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) + 
  xlab("Log(COVID-19 cases per capita)") + ylab("")    +  
  scale_fill_cyclical(values = c("#214D72", "#50BFC3")) + 
  scale_x_continuous(breaks = c(-12, -9, -6, -3, 0))  + coord_cartesian(xlim =c(-12, 0)) 

#coord_cartesian(xlim =c(0, -12)) 

library(gridExtra)
library(ggpubr)
library(grid)
library(egg)

jpeg("Figures/finalggridge_covdpercapita.jpg", units="in", width=14, height=7, res=300)

egg::ggarrange(f4, f3, ncol = 1, heights = c(0.89, 1))

dev.off()




##############################################################################
#        Main Figure 2(a) 
#############################################################################


jpeg("Figures/Top25Pop_LT100Pop.jpg", units="in", width=9, height=5.4, res=300)

r1 <- ggplot(full2[which(full2$Top25Pop_LT100Pop!="NA"),], 
             aes(x = SIP_Index_April1*100, 
                 y = Top25Pop_LT100Pop)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2, size=1, inherit.aes = TRUE) + 
  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.7,
    aes(fill= Top25Pop_LT100Pop)) +  
  xlab("% Stay in Place") + ylab("")    +  
  scale_x_continuous(breaks = c(20, 40, 60)) +
  scale_fill_manual(values =c("#214D72", "#50BFC3")) 
#+  coord_cartesian(xlim =c(0, 100)) 
dev.off()


jpeg("Figures/Top25_LT100_CovidCore.jpg", units="in", width=9, height=5.4, res=300)

r2 <- ggplot(full2[which(full2$Top25_LT100_CovidCore!="NA"),], 
             aes(x = SIP_Index_April1*100, 
                 y = Top25_LT100_CovidCore)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2, size=1, inherit.aes = TRUE) + 
  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.7,
    aes(fill= Top25_LT100_CovidCore)) +  
  xlab("% Stay in Place") + ylab("")    +  
  scale_x_continuous(breaks = c(20, 40, 60)) +
  scale_fill_manual(values =c("#214D72", "#50BFC3")) 
r2
dev.off()

jpeg("Figures/StayinPlaceRawComparison.jpg", units="in", width=9, height=7, res=300)

egg::ggarrange(r1, r2, ncol = 1, heights = c(0.89, 1))




#######################################################
# Figure 2 (b)                                       # 
#######################################################


travel.1 <- lm(SIP_Index_April1*100~ Top25Pop + stayathome_surveyweek + 
                 log_covidpercapitanew + pop.density + med.income +
                 perc.white + perc.less.hs + rep.perc16, data = full2)


travel.2 <- lm(SIP_Index_April1*100~ Top25CovidCore + stayathome_surveyweek + 
                 log_covidpercapitanew  + pop.density + med.income +
                 perc.white + perc.less.hs + rep.perc16, data = full2)


travel.3 <- lm(SIP_Index_April1*100~ logCovidCore + stayathome_surveyweek + 
                 log_covidpercapitanew  + pop.density + med.income +
                 perc.white + perc.less.hs + rep.perc16, data = full2)


m1 <- coef(summary(travel.1))
m2 <- coef(summary(travel.2))
m3 <- coef(summary(travel.3))

# Combined 
c3 <- numeric(length = 12)


c3[1] <- m1[9,1]*sd(full2$rep.perc16)
c3[2] <- m2[9,1]*sd(full2$rep.perc16)
c3[3] <- m3[9,1]*sd(full2$rep.perc16)

c3[4] <- m1[4,1]*sd(full2$log_covidpercapitanew, na.rm=T) 
c3[5] <- m2[4,1]*sd(full2$log_covidpercapitanew, na.rm=T) 
c3[6] <- m3[4,1]*sd(full2$log_covidpercapitanew, na.rm=T) 

c3[7] <- m1[3,1]*sd(full2$stayathome_surveyweek)
c3[8] <- m2[3,1]*sd(full2$stayathome_surveyweek)
c3[9] <- m3[3,1]*sd(full2$stayathome_surveyweek)

c3[10] <- m1[2,1]*sd(full2$Top25Pop, na.rm=T)
c3[11] <- m2[2,1]*sd(full2$Top25CovidCore, na.rm=T) 
c3[12] <- m3[2,1]*sd(full2$logCovidCore, na.rm=T) 



c3 <-as.data.frame((c3))
c3[1,c(2)] <- "#0F9FD6" 
c3[2,c(2)] <- "#307FE2"
c3[3,c(2)] <- "#1E3F63"

c3[4,c(2)] <- "#0F9FD6" 
c3[5,c(2)] <- "#307FE2"
c3[6,c(2)] <- "#1E3F63"

c3[7,c(2)]<- "#0F9FD6" 
c3[8,c(2)] <- "#307FE2"
c3[9,c(2)] <- "#1E3F63"

c3[10,c(2)] <- "#0F9FD6" 
c3[11,c(2)] <- "#307FE2"
c3[12,c(2)] <- "#1E3F63"

s3 <- numeric(length = 12)

s3[1] <- m1[9,2]*sd(full2$rep.perc16)
s3[2] <- m2[9,2]*sd(full2$rep.perc16)
s3[3] <- m3[9,2]*sd(full2$rep.perc16)

s3[4] <- m1[4,2]*sd(full2$log_covidpercapitanew, na.rm=T) 
s3[5] <- m2[4,2]*sd(full2$log_covidpercapitanew, na.rm=T) 
s3[6] <- m3[4,2]*sd(full2$log_covidpercapitanew, na.rm=T) 

s3[7] <- m1[3,2]*sd(full2$stayathome_surveyweek)
s3[8] <- m2[3,2]*sd(full2$stayathome_surveyweek)
s3[9] <- m3[3,2]*sd(full2$stayathome_surveyweek)

s3[10] <- m1[2,2]*sd(full2$Top25Pop,  na.rm=T)
s3[11] <- m2[2,2]*sd(full2$Top25CovidCore, na.rm=T) 
s3[12] <- m3[2,2]*sd(full2$logCovidCore, na.rm=T) 

yloc <- c(1.1,1.0,0.9, 2.1,2.0,1.9, 3.1,3.0, 2.9, 4.1,4.0, 3.9)

jpeg("Figures/SIP_CoefficientPlot.jpg", units="in", width=9, height=5, res=300)
p <- recordPlot()
plot.new()
plot(c3$`(c3)`, yloc, pch = 23, col=c3$V2, bg=c3$V2, 
     xlim = c(-5, 5), ylim = c(.5,4.5), xlab = "Impact of Standard Deviation Change", 
     ylab = "", yaxt = "n",axes=FALSE)
title("Variation in % Stay at Home")
abline(v = 0, col = "gray")
segments((c3$`(c3)` - (qnorm(0.95) * s3)), yloc, (c3$`(c3)` + (qnorm(0.95) * s3)), yloc, col = c3$V2, 
         lwd = 2)
segments((c3$`(c3)` - (qnorm(0.975) * s3)), yloc, (c3$`(c3)` + (qnorm(0.975) * s3)), yloc, col =  c3$V2, 
         lwd = 1)
grid()
axis(1)
text(-0.04,1.4,'% Trump Vote',pos=4)
text(-0.04,2.4,'log(COVID Cases per capita)', pos=4)
text(-0.04,3.4,'Stay Home Order',pos=4)
text(-0.04,4.4, 'Urban DMA Measure', pos=4)
text(2, .9, "Model 1: Top 25 DMA Indicator", col = "#0F9FD6" ,cex=.8,pos=4)
text(2, .7, "Model 2: Top 25 COVID Core", col = "#307FE2" ,cex=.8,pos=4)
text(2, .5, "Model 3: log(COVID Cases DMA) ", col = "#1E3F63",cex=.8,pos=4)
library(gridGraphics)
library(grid)
a<-    grid.grab()
dev.off() 

#######################################################
# Table 1                                             # 
#######################################################


stayhome.1a <- lm(COVIDBEHAVIOR_stayhome ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
stayhome.1b <- lm(COVIDBEHAVIOR_stayhome ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)

stayhome.2a <- lm(COVIDBEHAVIOR_stayhome ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
stayhome.2b <- lm(COVIDBEHAVIOR_stayhome ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)

maskoutside.1a <- lm(COVIDBEHAVIOR_maskoutside ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
maskoutside.1b <- lm(COVIDBEHAVIOR_maskoutside ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)

maskoutside.2a <- lm(COVIDBEHAVIOR_maskoutside ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
maskoutside.2b <- lm(COVIDBEHAVIOR_maskoutside ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)


table <- capture.output({ # Store the stargazer output in a string
  stargazer(stayhome.1a,stayhome.1b, stayhome.2a, stayhome.2b, maskoutside.1a, maskoutside.1b, maskoutside.2a, maskoutside.2b,
            dep.var.caption = "",
            dep.var.labels=c("Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("Top 25 DMA", "Top 25 Covid Core", "Stay at Home Order","log(COVID county cases per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            star.cutoffs = c(.1, .05,.01,.001),
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effects of Residing in Urban DMA in Self-Reported COVID-19 Social Distancing Behaviors",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#######################################################
#######################################################
# Appendix                                            # 
#######################################################
#######################################################



####  Figure S1


f3 <- ggplot(full2[which(full2$top25factor!="NA"),], 
             aes(x = log_dmacasespercapita, 
                 y = top25factor , 
                 fill = top25factor),  scale_fill_manual(values =c("#1D3752", "#50BFC3"))) +  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.3,
    aes(fill= top25factor)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) + 
  xlab("Log(COVID-19 cases in DMA County Core per capita)") + ylab("")    +  
  scale_fill_cyclical(values = c("#214D72", "#50BFC3")) + 
  scale_x_continuous(breaks = c(-12, -9, -6, -3, 0))  + coord_cartesian(xlim =c(-12, 0)) 


f4 <- ggplot(full2[which(full2$top25factor!="NA"),], 
             aes(x = log_covidpercapitanew, 
                 y = top25factor , 
                 fill = top25factor),  scale_fill_manual(values =c("#1D3752", "#50BFC3"))) +  
  theme_ridges() + theme(legend.position = "none", 
                         axis.title.x = element_text(color="black", size=14, face="bold"),
                         axis.text.y = element_text(color= "black", size=14, face="bold"),
                         axis.text.x = element_text(color= "black", size=14, face="bold"),
                         axis.title.y = element_text(color="black", size=14, face="bold")) + 
  geom_density_ridges(
    jittered_points = TRUE, 
    position = position_points_jitter(width = 0.05, height = 0),
    point_shape = "|", point_size = 3,  point_alpha = 1, alpha = 0.7,
    aes(fill= top25factor)) +  
  stat_density_ridges(quantile_lines = TRUE, quantiles = 2) + 
  xlab("Log(COVID-19 cases per capita)") + ylab("")    +  
  scale_fill_cyclical(values = c("#214D72", "#50BFC3")) + 
  scale_x_continuous(breaks = c(-12, -9, -6, -3, 0))  + coord_cartesian(xlim =c(-12, 0)) 


egg::ggarrange(f4, f3, ncol = 1, heights = c(0.89, 1))


#### Table S2

t1<- t.test(ind$COVIDBEHAVIOR_maskoutside ~   ind$Top25Pop) 
t1
summary(ind$COVIDBEHAVIOR_maskoutside)
t2 <-t.test(ind$COVIDBEHAVIOR_stayhome ~    ind$Top25Pop) 
t2
summary(ind$COVIDBEHAVIOR_stayhome)
t3<- t.test(ind$COVIDCONCERNR1  ~ind$Top25Pop) 
t3
summary(ind$COVIDCONCERNR1)
t4 <-t.test(ind$Dem ~      ind$Top25Pop) 
t4
summary(ind$Dem)
t5<- t.test(ind$Rep~     ind$Top25Pop) 
t5
summary(ind$Rep)
t6<- t.test(ind$Age~  ind$Top25Pop) 
t6
summary(ind$Age)
t7<- t.test(ind$ChurchWeekly~  ind$Top25Pop) 
t7
summary(ind$ChurchWeekly)
t8<- t.test(ind$CHILD~  ind$Top25Pop) 
t8
summary(ind$CHILD)
t9<- t.test(ind$OLDPARENTS~  ind$Top25Pop) 
t9
summary(ind$OLDPARENTS)
t10<- t.test(ind$HSEduc~  ind$Top25Pop) 
t10
summary(ind$HSEduc)
t11<- t.test(ind$CollegeEduc~  ind$Top25Pop) 
t11
summary(ind$CollegeEduc)
t12<- t.test(ind$Male~  ind$Top25Pop) 
t12
summary(ind$Male)
t13
t13<- t.test(ind$UNEMPLOYED~  ind$Top25Pop) 
summary(ind$UNEMPLOYED)

#### Table S3 

sort(full2$Population,decreasing=TRUE,index.return=TRUE)
sort(full2$DMACoreCases,decreasing=TRUE,index.return=TRUE)


#### Table S4 



travel.1 <- lm(SIP_Index_April1*100~ Top25Pop + stayathome_surveyweek + 
                 log_covidpercapitanew + pop.density + med.income +
                 perc.white + perc.less.hs + rep.perc16, data = full2)
summary(travel.1)


travel.1b <- lm(YOYChange_Apr1*100~ Top25Pop + stayathome_surveyweek + 
                  log_covidpercapitanew + pop.density + med.income +
                  perc.white + perc.less.hs  + rep.perc16, data = full2)
summary(travel.1b)


travel.2 <- lm(SIP_Index_April1*100~ Top25CovidCore + stayathome_surveyweek + 
                 log_covidpercapitanew  + pop.density + med.income +
                 perc.white + perc.less.hs + rep.perc16, data = full2)
summary(travel.2)

travel.2b <- lm(YOYChange_Apr1*100~ Top25CovidCore + stayathome_surveyweek + 
                  log_covidpercapitanew  + pop.density + med.income +
                  perc.white + perc.less.hs + rep.perc16, data = full2)
summary(travel.2b)


travel.3 <- lm(SIP_Index_April1*100~ logCovidCore + stayathome_surveyweek + 
                 log_covidpercapitanew  + pop.density + med.income +
                 perc.white + perc.less.hs + rep.perc16, data = full2)
summary(travel.3)

travel.3b <- lm(YOYChange_Apr1*100~ logCovidCore + stayathome_surveyweek + 
                  log_covidpercapitanew  + pop.density + med.income +
                  perc.white + perc.less.hs + rep.perc16 , data = full2)
summary(travel.3b)



table <- capture.output({stargazer(travel.1, travel.1b, travel.2, travel.2b, travel.3, travel.3b,
                                   dep.var.labels=c("% Stay at Home", "Change in % Stay at Home", "% Stay at Home", "Change in % Stay at Home", "% Stay at Home", "Change in % Stay at Home"),
                                   covariate.labels = c("Top 25 DMA","Top 25 COVID Core","log(Cases in DMA County Core)","Stay at Home Order","log(COVID cases in County per capita)","Pop. Density","Median Income","% White","% HS Educ. or less","% Trump Vote"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   no.space=TRUE,
                                   star.char = c("+", "*", "**", "***"), 
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   star.cutoffs = c(.1, .05,.01,.001),
                                   notes.append = F,
                                   notes.align="l",
                                   digits=1,
                                   align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#### Table S5 

stravel.1 <- lm(SIP_Index_April1*100~ Top25Pop + 
                  log_covidpercapitanew  + pop.density + med.income +
                  perc.white + perc.less.hs + rep.perc16 + as.factor(statefp), data = full2)
summary(stravel.1)

stravel.1b <- lm(YOYChange_Apr1*100~ Top25Pop + 
                   log_covidpercapitanew  + pop.density + med.income +
                   perc.white + perc.less.hs  + rep.perc16 + as.factor(statefp), data = full2)
summary(stravel.1b)

stravel.2 <- lm(SIP_Index_April1*100~ Top25CovidCore + 
                  log_covidpercapitanew  + pop.density + med.income +
                  perc.white + perc.less.hs + rep.perc16 + as.factor(statefp), data = full2)
summary(stravel.2)

stravel.2b <- lm(YOYChange_Apr1*100~ Top25CovidCore + 
                   log_covidpercapitanew  + pop.density + med.income +
                   perc.white + perc.less.hs + rep.perc16 + as.factor(statefp), data = full2)
summary(stravel.2b)

stravel.3 <- lm(SIP_Index_April1*100~ logCovidCore + 
                  log_covidpercapitanew  + pop.density + med.income +
                  perc.white + perc.less.hs + rep.perc16 + as.factor(statefp), data = full2)
summary(stravel.3)

stravel.3b <- lm(YOYChange_Apr1*100~ logCovidCore  + 
                   log_covidpercapitanew  + pop.density + med.income +
                   perc.white + perc.less.hs + rep.perc16 + as.factor(statefp), data = full2)
summary(stravel.3b)



table <- capture.output({stargazer(stravel.1, stravel.1b, stravel.2, stravel.2b, stravel.3, stravel.3b,
                                   dep.var.labels=c("% Stay at Home", "Change in % Stay at Home", "% Stay at Home", "Change in % Stay at Home"),
                                   covariate.labels = c("Top 25 DMA", "Top 25 Covid Core", "log(Cases in DMA County Core)","log(COVID cases in County per capita)","Pop. Density","Median Income","% White","% HS Educ. or less","% Trump Vote"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = c("statefp"),
                                   omit.labels = c("State Fixed Effects"),
                                   no.space=TRUE,
                                   star.char = c("+", "*", "**", "***"), 
                                   star.cutoffs = c(.1, .05,.01,.001),
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=1,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#### Table S6


library(cobalt)
library(backports)

t1<- t.test(full$pop.density ~     full$Top25Pop.x) 
t3 <-t.test(full$perc.white ~    full$Top25Pop.x) 
t4<- t.test(full$rep.perc16.x  ~      full$Top25Pop.x) 
t2 <-t.test(full$med.income ~      full$Top25Pop.x) 
t5<- t.test(full$perc.less.hs~     full$Top25Pop.x) 
t6<- t.test(full$log_covidpercapitanew~  full$Top25Pop.x) 
t7 <- t.test(full$perc.unemployed ~ full$Top25Pop.x)

# Balance table 

dtm <- full[!is.na(full$Top25Pop),]
covs <- subset(dtm, select = c(Top25Pop, pop.density, perc.white, rep.perc16.x, med.income, perc.less.hs, log_covidpercapita))

bal.tab(Top25Pop ~ covs, data = dtm,  
        disp = c("means"), un = TRUE) 

pval <- c(t1[3], t2[3] ,t3[3],t4[3],t5[3],t6[3])

#### Table S7

#full$Top25Pop <- full$Top25Pop.x
#full$LT100Pop <- full$LT100Pop.x
#full$Top25CovidCore <- full$Top25CovidCore.x
#full$LT100CovidCore <- full$LT100CovidCore.x
#full$DMACoreCases <- full$DMACoreCases.x 


m1 <- lm(SIP_Index_April1*100~ Top25Pop + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full)

m1s <- lm(SIP_Index_April1*100~ Top25Pop + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full)

m2 <- lm(YOYChange_Apr1*100~ Top25Pop + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full)

m2s <- lm(YOYChange_Apr1*100~ Top25Pop + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full)


m3 <- lm(SIP_Index_April1*100~ Top25Pop + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek, data=full[full$Top25Pop==1 | full$LT100Pop==1,]) 

m3s <- lm(SIP_Index_April1*100~ Top25Pop + log_covidpercapitanew  + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp), data=full[full$Top25Pop==1 | full$LT100Pop==1,])

m4 <- lm(YOYChange_Apr1*100~ Top25Pop + log_covidpercapitanew  + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek, data=full[full$Top25Pop==1 | full$LT100Pop==1,])

m4s <- lm(YOYChange_Apr1*100~ Top25Pop + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full[full$Top25Pop==1 | full$LT100Pop==1,]) 


table <- capture.output({stargazer(m1, m1s, m2, m2s, m3, m3s, m4, m4s, 
                                   dep.var.labels=c("% Stay at Home", "% Stay at Home", "Change in % Stay at Home", "Change in % Stay at Home", "% Stay at Home", "% Stay at Home", "Change in % Stay at Home", "Change in % Stay at Home"),
                                   covariate.labels = c("Top 25 DMA", "log(COVID cases in County per capita)", "Pop. Density","Median Income","% White", "% HS Educ. or less","% Trump Vote", "Stay at Home Order"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   omit = c("statefp"),
                                   omit.labels = c("State Fixed Effects"),
                                   no.space=TRUE,
                                   star.cutoffs = c(.1, .05,.01,.001),
                                   star.char = c("+", "*", "**", "***"), 
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=1,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#### Table S8 


m1 <- lm(SIP_Index_April1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full) 

m1s <- lm(SIP_Index_April1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full) 

m2 <- lm(YOYChange_Apr1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full)

m2s <- lm(YOYChange_Apr1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full)

m3 <- lm(SIP_Index_April1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full[full$Top25CovidCore==1 | full$LT100CovidCore==1,])

m3s <-lm(SIP_Index_April1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full[full$Top25CovidCore==1 | full$LT100CovidCore==1,]) 

m4 <-lm(YOYChange_Apr1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
          perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full[full$Top25CovidCore==1 | full$LT100CovidCore==1,])  

m4s <- lm(YOYChange_Apr1*100~ Top25CovidCore + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full[full$Top25CovidCore==1 | full$LT100CovidCore==1,]) 

table <- capture.output({stargazer(m1, m1s, m2, m2s, m3, m3s, m4, m4s, 
                                   dep.var.labels=c("% Stay at Home", "Change in % Stay at Home",  "% Stay at Home",  "Change in % Stay at Home"),
                                   covariate.labels = c("Top 25 Covid Core", "log(COVID cases in County per capita)", "Pop. Density","Median Income","% White", "% HS Educ. or less","% Trump Vote", "Stay at Home Order"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   no.space=TRUE,
                                   omit = c("statefp"),
                                   omit.labels = c("State Fixed Effects"),
                                   star.cutoffs = c(.1, .05,.01,.001),
                                   star.char = c("+", "*", "**", "***"), 
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=1,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#### Table S9 

m1 <- lm(SIP_Index_April1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full)

m1s <- lm(SIP_Index_April1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full)

m2 <- lm(YOYChange_Apr1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full)

m2s <- lm(YOYChange_Apr1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full)


m3 <- lm(SIP_Index_April1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full[full$Top25Pop==1 | full$LT100Pop==1,])

m3s <- lm(SIP_Index_April1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full[full$Top25Pop==1 | full$LT100Pop==1,])

m4 <- lm(YOYChange_Apr1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
           perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek,data=full[full$Top25Pop==1 | full$LT100Pop==1,])

m4s <- lm(YOYChange_Apr1*100~ log(DMACoreCases) + log_covidpercapitanew + pop.density + med.income +
            perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp),data=full[full$Top25Pop==1 | full$LT100Pop==1,]) 

table <- capture.output({stargazer(m1, m1s, m2, m2s, m3, m3s, m4, m4s, 
                                   dep.var.labels=c("% Stay at Home", "Change in % Stay at Home",  "% Stay at Home",  "Change in % Stay at Home"),
                                   covariate.labels = c("log(DMA Core Cases)", "log(COVID cases in County per capita)", "Pop. Density","Median Income","% White", "% HS Educ. or less","% Trump Vote", "Stay at Home Order"),
                                   omit = c("statefp"),
                                   omit.labels = c("State Fixed Effects"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   no.space=TRUE,
                                   star.cutoffs = c(.1, .05,.01,.001),
                                   star.char = c("+", "*", "**", "***"), 
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=1,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


#### Table S11 : make sure to download county_adjacency.csv file. 

# extracting the list of county fips codes for Top25 (most popuplous) DMA 

fipslist <- (full$fips[full$Top25Pop==1])
fipslist <- as.data.frame(full$fips[full$Top25Pop==1])
colnames(fipslist) <- "fips"
fipslist <-fipslist[!is.na(fipslist$fips),]
fipslist <- as.data.frame(fipslist)
adj <- read.csv("county_adjacency.csv")
colnames(fipslist) <- "fips"


adj_wide <- dcast(adj, county ~ neighboring.county, value.var = "neighboring.county")


colnames(adj_wide)[1] <- "fips" 

adj_top25 <- merge(fipslist, adj_wide, all.x=TRUE, by="fips")


# transpose table 

adj_top25_t <- t(adj_top25)
adj_top25_t <- as.data.frame(adj_top25_t)
colnames(adj_top25_t) <- adj_top25_t[1,]
adj_top25_t <- adj_top25_t[-1, ]

# extract all variable names when the sum of each row is greater than 1 

adj_top25_t$rowtotal <- rowSums(adj_top25_t)  
adj_top25_final <- adj_top25_t[which(adj_top25_t$rowtotal > 0), ]

# extracting row numbers (that would give us the list of fips code that are adjacent counties)

r <- rownames(adj_top25_final)
r <- as.matrix(r)
r <- as.data.frame(r)
colnames(r)[1] <- "fips" 
r$top25adj_dummy <- 1 
r$fips <- as.numeric(as.character(r$fips))

#write.csv(r, "Data/top25adj_list.csv" )

# now merge with full 

#r <- read.csv("Data/top25adj_list.csv")

dtm <- merge(full, r, all= TRUE, by = "fips")
#dtm$top25adj_dummy[is.na(dtm$top25adj_dummy)] <- 0 

# generate new treatmet variable 
# now existing top25pop is the treatment group (1) and
# adjacenet counties of top 25 counties are control (0)
# there isa total of 387 treated counties and 202 adjacenet counties 
# HOWEVER, there are overlaps between these counties
# 116 counties overlap. we exclude these counties from treatment group. 
# that gives us 271 treated counties and 202 adjacement counties. 

dtm$top25adj_dummy[is.na(dtm$top25adj_dummy)] <- 0 
# exlcude any adjacent FIPs county code that belongs to Top25Pop
dtm$top25adj_dummy[dtm$top25adj_dummy==1 & dtm$Top25Pop==1] <- 0 


dtm$top25_t_adj_c  <- NA 
dtm$top25_t_adj_c[dtm$Top25Pop==1] <- 1 
dtm$top25_t_adj_c[dtm$top25adj_dummy==1] <- 0 

dtm_adj <- dtm[!is.na(dtm$top25_t_adj_c),]

adj1  <- lm(SIP_Index_April1*100~ top25_t_adj_c +   
              log_covidpercapitanew + pop.density + med.income +
              perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek, data = dtm_adj)
summary(adj1)

#   Add State FE

adj1s <- lm(SIP_Index_April1*100~ top25_t_adj_c +  
              log_covidpercapitanew + pop.density + med.income +
              perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp), data = dtm)
summary(adj1s )


# adjacency balance table 

dtm_adj$top25_t_adj_c <- as.factor( dtm_adj$top25_t_adj_c)
t1<- t.test(dtm_adj$pop.density ~ dtm_adj$top25_t_adj_c) 
t2 <- t.test(dtm_adj$perc.white ~ dtm_adj$top25_t_adj_c) 
t3 <-t.test(dtm_adj$rep.perc16.x ~ dtm_adj$top25_t_adj_c) 
t4<- t.test(dtm_adj$med.income ~ dtm_adj$top25_t_adj_c) 
t5<- t.test(dtm_adj$perc.less.hs ~ dtm_adj$top25_t_adj_c) 
t6 <- t.test(dtm_adj$log_covidpercapitanew ~  dtm_adj$top25_t_adj_c) 
t7 <- t.test(dtm_adj$perc.unemployed~  dtm_adj$top25_t_adj_c) 

# Balance table 
covs <- subset(dtm_adj, select = c(pop.density, perc.white, rep.perc16.x, med.income, perc.less.hs, perc.unemployed, log_covidpercapitanew))

bal.tab(top25_t_adj_c ~ covs, data = dtm_adj, 
        disp = c("means"), un = FALSE, 
        stats = c("mean.diffs"))

pval <- c(t1[3], t2[3] ,t3[3],t4[3],t5[3],t7[3], t6[3])
pval


#     Adjacency (Analysis for All Rural Counties Adjacent to a rural county in a Top COVID 25)

# extracting the list of county fips codes for Top25 (most COVID impacted) DMA 
# total of 387 rural counties belong to Top 25 DMA


fipslist <- (full$fips[full$Top25CovidCore==1])
fipslist <- as.data.frame(full$fips[full$Top25CovidCore==1])
colnames(fipslist) <- "fips"
fipslist <-fipslist[!is.na(fipslist$fips),]
fipslist <- as.data.frame(fipslist)
adj <- read.csv("Data/county_adjacency.csv")
colnames(fipslist) <- "fips"



adj_wide <- dcast(adj, county ~ neighboring.county, value.var = "neighboring.county")

colnames(fipslist) <- "fips"
colnames(adj_wide)[1] <- "fips" 

adj_top25 <- merge(fipslist, adj_wide, all.x=TRUE, by="fips")


# transpose table 

adj_top25_t <- t(adj_top25)
adj_top25_t <- as.data.frame(adj_top25_t)
colnames(adj_top25_t) <- adj_top25_t[1,]
adj_top25_t <- adj_top25_t[-1, ]

# extract all variable names when the sum of each row is greater than 1 

adj_top25_t$rowtotal <- rowSums(adj_top25_t)  
adj_top25_final <- adj_top25_t[which(adj_top25_t$rowtotal > 0), ]

# extracting row numbers (that would give us the list of fips code that are adjacent counties)

r <- rownames(adj_top25_final)
r <- as.data.frame(r)
colnames(r) <- "fips" 
r$top25adj_dummy <- 1 


# now merge with full 

dtm <- merge(full, r, all.x= TRUE, by = "fips")
dtm$top25adj_dummy[is.na(dtm$top25adj_dummy)] <- 0 
dtm$top25adj_dummy[dtm$top25adj_dummy==1 & dtm$Top25CovidCore==1] <- 0 

# generate new treatmet variable 
# now existing top25pop is the treatment group (1) and
# adjacent counties of top 25 counties are control (0)

dtm$top25_t_adj_c  <- NA 
dtm$top25_t_adj_c[dtm$Top25CovidCore==1] <- 1 
dtm$top25_t_adj_c[dtm$top25adj_dummy==1] <- 0 
# this gives us a total of 374 treated counties and 68 adjacenet counties 

dtm_adj <- dtm[!is.na(dtm$top25_t_adj_c),]


adj2  <- lm(SIP_Index_April1*100~ top25_t_adj_c +   
              log_covidpercapitanew + pop.density + med.income +
              perc.white + perc.less.hs + rep.perc16.x + stayathome_surveyweek, data = dtm_adj)
summary(adj2)

#   Add State FE

adj2s <- lm(SIP_Index_April1*100~ top25_t_adj_c +  
              log_covidpercapitanew + pop.density + med.income +
              perc.white + perc.less.hs + rep.perc16.x + as.factor(statefp), data = dtm_adj)
summary(adj2s)

table <- capture.output({stargazer(adj1, adj1s, adj2, adj2s, 
                                   dep.var.labels=c("% Stay at Home"),
                                   covariate.labels = c("Top 25 Rural Counties vs Adjacent Rural Counties", "log(COVID cases in County per capita)", "Pop. Density","Median Income","% White", "% HS Educ. or less","% Trump Vote", "Stay at Home Order"),
                                   omit = c("statefp"),
                                   omit.labels = c("State Fixed Effects"),
                                   omit.stat=c("adj.rsq","LL","ser","f"),
                                   no.space=TRUE,
                                   star.char = c("+", "*", "**", "***"), 
                                   notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
                                   notes.append = F,
                                   notes.align="l",
                                   digits=1,
                                   align = TRUE) 
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


####  Table S12, S13, S14


MediaUsage.2 <- lm(NEWSHR ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
MediaUsage.3 <- lm(NEWSHR ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
MediaUsage.3nl <- lm(NEWSHR ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

MediaUsage.4 <- lm(NEWSHR ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
MediaUsage.5 <- lm(NEWSHR ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
MediaUsage.5nl <- lm(NEWSHR ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

# Prefer to Watch More local news?

ind$localnewswatch <- NA
ind$localnewswatch[ind$MEDIA2==2] <- 1
ind$localnewswatch[ind$MEDIA2==3] <- 1
ind$localnewswatch[ind$MEDIA2==1] <- 0
ind$localnewswatch[ind$MEDIA2==4] <- 0
datlocal$localnewswatch <- NA
datlocal$localnewswatch[datlocal$MEDIA2==2] <- 1
datlocal$localnewswatch[datlocal$MEDIA2==3] <- 1
datlocal$localnewswatch[datlocal$MEDIA2==1] <- 0
datlocal$localnewswatch[datlocal$MEDIA2==4] <- 0

datnolocal$localnewswatch <- NA
datnolocal$localnewswatch[datnolocal$MEDIA2==2] <- 1
datnolocal$localnewswatch[datnolocal$MEDIA2==3] <- 1
datnolocal$localnewswatch[datnolocal$MEDIA2==1] <- 0
datnolocal$localnewswatch[datnolocal$MEDIA2==4] <- 0

MoreLocal.2 <- lm(localnewswatch ~ Top25Pop   + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
MoreLocal.3 <- lm(localnewswatch ~ Top25Pop   + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
MoreLocal.3nl <- lm(localnewswatch ~ Top25Pop   + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

MoreLocal.4 <- lm(localnewswatch ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
MoreLocal.5 <- lm(localnewswatch ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
MoreLocal.5nl <- lm(localnewswatch ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

# Approve National News?

NatNewsApp.2 <- lm(COVIDSAT_nationalnews ~ Top25Pop  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
NatNewsApp.3 <- lm(COVIDSAT_nationalnews ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
NatNewsApp.3nl <- lm(COVIDSAT_nationalnews ~ Top25Pop  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

NatNewsApp.4 <- lm(COVIDSAT_nationalnews ~   log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
NatNewsApp.5 <- lm(COVIDSAT_nationalnews ~   log(DMACoreCases)+ Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
NatNewsApp.5nl <- lm(COVIDSAT_nationalnews ~ log(DMACoreCases)+ Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

# Approve Local News (top 25 = less likley to approve )

LocalNewsApp.2 <- lm(COVIDSAT_localnews ~ Top25Pop  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
LocalNewsApp.3 <- lm(COVIDSAT_localnews ~ Top25Pop  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
LocalNewsApp.3nl <- lm(COVIDSAT_localnews ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

LocalNewsApp.4 <- lm(COVIDSAT_localnews ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
LocalNewsApp.5 <- lm(COVIDSAT_localnews ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
LocalNewsApp.5nl <- lm(COVIDSAT_localnews ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

# total number of local media - No difference! 

LocalNewsTotal.2 <- lm(LOCALNEWS_sourcetotal ~   Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
LocalNewsTotal.3 <- lm(LOCALNEWS_sourcetotal ~   Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
LocalNewsTotal.3nl <- lm(LOCALNEWS_sourcetotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

LocalNewsTotal.4 <- lm(LOCALNEWS_sourcetotal ~   log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
LocalNewsTotal.5 <- lm(LOCALNEWS_sourcetotal ~   log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
LocalNewsTotal.5nl <- lm(LOCALNEWS_sourcetotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

# Fox News Consumption - no diff 

Fox.2 <- lm(FOXtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
Fox.3 <- lm(FOXtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
Fox.3nl <- lm(FOXtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

Fox.4 <- lm(FOXtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
Fox.5 <- lm(FOXtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
Fox.5nl <- lm(FOXtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

#MSNBC Consumption  - no diff 
MSNBC.2 <- lm(MSNBCtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
MSNBC.3 <- lm(MSNBCtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
MSNBC.3nl <- lm(MSNBCtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

MSNBC.4 <- lm(MSNBCtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
MSNBC.5 <- lm(MSNBCtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
MSNBC.5nl <- lm(MSNBCtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)


#CNN  Consumption  - no diff 
CNN.2 <- lm(CNNtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
CNN.3 <- lm(CNNtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
CNN.3nl <- lm(CNNtotal ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

CNN.4 <- lm(CNNtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
CNN.5 <- lm(CNNtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
CNN.5nl <- lm(CNNtotal ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

#network 
NetworkEve.2 <- lm(NetworkEvening ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
NetworkEve.3 <- lm(NetworkEvening ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
NetworkEve.3nl <- lm(NetworkEvening ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)

NetworkEve.4 <- lm(NetworkEvening ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
NetworkEve.5 <- lm(NetworkEvening ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datlocal)
NetworkEve.5nl <- lm(NetworkEvening ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=datnolocal)


##################################################
#     TABLE S12 
#     Full Sample
##################################################

# full manipluation checks for appendix 
table <- capture.output({ # Store the stargazer output in a string
  stargazer(MediaUsage.2,MoreLocal.2,NatNewsApp.2,LocalNewsApp.2,LocalNewsTotal.2, Fox.2, MSNBC.2, CNN.2, NetworkEve.2,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Attitudes Toward Media and Media Consumption Patterns by Media Market Type",
            dep.var.labels=c("Media Usage (Hrs)","Watch Local News?","National News App.","Local News App.","Num. Local News Sources", "Num. of Fox News Programs", "Num. of MSNBC Programs", "Num. of CNN Programs", "Num. of Evening News Programs"),
            covariate.labels = c("Top 25 DMA","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per capita)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.cutoffs = c(.1, .05,.01,.001),
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

##################################################
#     TABLE S13
# full manipluation checks for appendix - only with local news consumers 
##################################################

table <- capture.output({ # Store the stargazer output in a string
  stargazer(MediaUsage.3,NatNewsApp.3,LocalNewsApp.3,LocalNewsTotal.3, Fox.3, MSNBC.3, CNN.3, NetworkEve.3,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Attitudes Toward Media and Media Consumption Patterns by Media Market Type: Local News Watchers",
            dep.var.labels=c("Media Usage (Hrs)","National News App.","Local News App.","Num of Local News Sources", "Num of Fox News", "Num of MSNBC", "Num of CNN", "Num of Evening News"),
            covariate.labels = c("Top 25 DMA","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per capita)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.char = c("+", "*", "**", "***"), 
            star.cutoffs = c(.1, .05,.01,.001),
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

##################################################
#     TABLE S14
#     only with non-local news consumers 
##################################################

table <- capture.output({ # Store the stargazer output in a string
  stargazer(MediaUsage.3nl,NatNewsApp.3nl,LocalNewsApp.3nl,LocalNewsTotal.3nl, Fox.3nl, MSNBC.3nl, CNN.3nl, NetworkEve.3nl,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Attitudes Toward Media and Media Consumption Patterns by Media Market Type: Non-Local Watchers",
            dep.var.labels=c("Media Usage (Hrs)","National News App.","Local News App.","Num of Local News Sources", "Num of Fox News", "Num of MSNBC", "Num of CNN", "Num of Evening News"),
            covariate.labels = c("Top 25 DMA","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per capita)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.cutoffs = c(.1, .05,.01,.001),
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

##################################################
#     TABLE S15
#     Full Sample 
#     using log(DMA cases) as a proxy 
##################################################

table <- capture.output({ # Store the stargazer output in a string
  stargazer(MediaUsage.4,MoreLocal.4,NatNewsApp.4,LocalNewsApp.4,LocalNewsTotal.4, Fox.4, MSNBC.4, CNN.4, NetworkEve.4,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Attitudes Toward Media and Media Consumption Patterns by COVID cases in DMA: All",
            dep.var.labels=c("Media Usage (Hrs)","Watch Local News?","National News App.","Local News App.","Num of Local News Sources", "Num of Fox News", "Num of MSNBC", "Num of CNN", "Num of Evening News"),
            covariate.labels = c("log(DMA Cty Cases)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per capita)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.cutoffs = c(.1, .05,.01,.001),
            
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

##################################################
#     TABLE S16
#     Local News Watchers Sample 
#     using log(DMA cases) as a proxy 
##################################################

table <- capture.output({ # Store the stargazer output in a string
  stargazer(MediaUsage.5,NatNewsApp.5,LocalNewsApp.5,LocalNewsTotal.5, Fox.5, MSNBC.5, CNN.5, NetworkEve.5,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Attitudes Toward Media and Media Consumption Patterns by COVID Cases in DMA: Local News Watchers",
            dep.var.labels=c("Media Usage (Hrs)","National News App.","Local News App.","Num of Local News Sources", "Num of Fox News", "Num of MSNBC", "Num of CNN", "Num of Evening News"),
            covariate.labels = c("log(DMA Cty Cases)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per capita)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.char = c("+", "*", "**", "***"), 
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

##################################################
#     TABLE S17
#     No Local News Watchers Sample 
#     using log(DMA cases) as a proxy 
##################################################

table <- capture.output({ # Store the stargazer output in a string
  stargazer(MediaUsage.5nl,NatNewsApp.5nl,LocalNewsApp.5nl,LocalNewsTotal.5nl, Fox.5nl, MSNBC.5nl, CNN.5nl, NetworkEve.5nl,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Attitudes Toward Media and Media Consumption Patterns by COVID Cases in DMA: Non-Local Watchers",
            dep.var.labels=c("Media Usage (Hrs)","National News App.","Local News App.","Num of Local News Sources", "Num of Fox News", "Num of MSNBC", "Num of CNN", "Num of Evening News"),
            covariate.labels = c("log(DMA Cty Cases)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per captia)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.cutoffs = c(.1, .05,.01,.001),
            
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


##################################################
#     TABLE S18
#     Comparing Local News Watchers    
##################################################


ind$localnewswatch <- NA
ind$localnewswatch[ind$MEDIA2==2] <- 1
ind$localnewswatch[ind$MEDIA2==3] <- 1
ind$localnewswatch[ind$MEDIA2==1] <- 0
ind$localnewswatch[ind$MEDIA2==4] <- 0

WatchLocal.1 <- lm(localnewswatch ~ Top25Pop + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)
WatchLocal.2 <- lm(localnewswatch ~ log(DMACoreCases) + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male +  UNEMPLOYED + log_covidpercapitanew,data=ind)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(WatchLocal.1,WatchLocal.2,
            dep.var.caption = "",
            digits=2,
            digits.extra = 1,
            title="Comparing Local News Watchers between Groups",
            covariate.labels = c("Top 25 DMA","log(DMA Cty Cases)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male", "Unemployed", "log(COVID cases in County per capita)"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.cutoffs = c(.1, .05,.01,.001),
            
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


#### Table S19: TABLE 1 Using Logistic Regression 

stayhome.1a <- glm(COVIDBEHAVIOR_stayhome ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal, family = "binomial")
stayhome.1b <- glm(COVIDBEHAVIOR_stayhome ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal, family = "binomial")

stayhome.2a <- glm(COVIDBEHAVIOR_stayhome ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal, family = "binomial")
stayhome.2b <- glm(COVIDBEHAVIOR_stayhome ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal, family = "binomial")

maskoutside.1a <- glm(COVIDBEHAVIOR_maskoutside ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal, family = "binomial")
maskoutside.1b <- glm(COVIDBEHAVIOR_maskoutside ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal, family = "binomial")

maskoutside.2a <- glm(COVIDBEHAVIOR_maskoutside ~ Top25Pop + stayathome_surveyweek + log_covidpercapitanew  + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal, family = "binomial")
maskoutside.2b <- glm(COVIDBEHAVIOR_maskoutside ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal, family = "binomial")


table <- capture.output({ # Store the stargazer output in a string
  stargazer(stayhome.1a,stayhome.1b, stayhome.2a, stayhome.2b, maskoutside.1a, maskoutside.1b, maskoutside.2a, maskoutside.2b,
            dep.var.caption = "",
            dep.var.labels=c("Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("Top 25 DMA", "Top 25 Covid Core", "Stay at Home Order","log(COVID county cases per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local", "No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effects of Residing in Urban DMA in Self-Reported COVID-19 Social Distancing Behaviors",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


####   Table S20: Using DMACoreCases 

stayhome.1c <- lm(COVIDBEHAVIOR_stayhome ~ log(DMACoreCases) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
stayhome.2c <- lm(COVIDBEHAVIOR_stayhome ~ log(DMACoreCases) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
maskoutside.1c <- lm(COVIDBEHAVIOR_maskoutside ~ log(DMACoreCases) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
maskoutside.2c <- lm(COVIDBEHAVIOR_maskoutside ~ log(DMACoreCases) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(stayhome.1c, stayhome.2c,  maskoutside.1c, maskoutside.2c,
            dep.var.caption = "",
            dep.var.labels=c("Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("log(DMA Core Cases)", "Stay at Home Order","log(COVID county cases per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effects of Residing in Urban DMA in Self-Reported COVID-19 Social Distancing Behaviors",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

####  TABLE S21  Using DMACoreCAses per capita 

datlocal$DMACoreCasesPerCapita <- datlocal$DMACoreCases / datlocal$DMACorePopulation
datnolocal$DMACoreCasesPerCapita <- datnolocal$DMACoreCases / datnolocal$DMACorePopulation

stayhome.1c <- lm(COVIDBEHAVIOR_stayhome ~ log(DMACoreCasesPerCapita) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
stayhome.2c <- lm(COVIDBEHAVIOR_stayhome ~ log(DMACoreCasesPerCapita) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
maskoutside.1c <- lm(COVIDBEHAVIOR_maskoutside ~ log(DMACoreCasesPerCapita) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
maskoutside.2c <- lm(COVIDBEHAVIOR_maskoutside ~ log(DMACoreCasesPerCapita) + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(stayhome.1c, stayhome.2c,  maskoutside.1c, maskoutside.2c,
            dep.var.caption = "",
            dep.var.labels=c("Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("log(DMA Core Cases per capita)", "Stay at Home Order","log(COVID county cases per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effects of Residing in Urban DMA in Self-Reported COVID-19 Social Distancing Behaviors",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


#### Table S22: Using DMA death

stayhome.1c <- lm(COVIDBEHAVIOR_stayhome ~ DMACoreDeaths + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
stayhome.2c <- lm(COVIDBEHAVIOR_stayhome ~ DMACoreDeaths + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
maskoutside.1c <- lm(COVIDBEHAVIOR_maskoutside ~ DMACoreDeaths + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
maskoutside.2c <- lm(COVIDBEHAVIOR_maskoutside ~ DMACoreDeaths + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(stayhome.1c, stayhome.2c,  maskoutside.1c, maskoutside.2c,
            dep.var.caption = "",
            dep.var.labels=c("Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("DMA Core COVID Deaths", "Stay at Home Order","log(COVID county cases per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effects of Residing in Urban DMA in Self-Reported COVID-19 Social Distancing Behaviors",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


####  Table S23: Using DMA Core death per capita 

datlocal$DMACoreDeathPerCapita <- datlocal$DMACoreDeaths / datlocal$DMACorePopulation
datnolocal$DMACoreDeathPerCapita <- datnolocal$DMACoreDeaths / datnolocal$DMACorePopulation

stayhome.1c <- lm(COVIDBEHAVIOR_stayhome ~ DMACoreDeathPerCapita + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
stayhome.2c <- lm(COVIDBEHAVIOR_stayhome ~ DMACoreDeathPerCapita + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
maskoutside.1c <- lm(COVIDBEHAVIOR_maskoutside ~ DMACoreDeathPerCapita + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
maskoutside.2c <- lm(COVIDBEHAVIOR_maskoutside ~ DMACoreDeathPerCapita + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(stayhome.1c, stayhome.2c,  maskoutside.1c, maskoutside.2c,
            dep.var.caption = "",
            dep.var.labels=c("Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("DMA Core COVID Death Per Capita", "Stay at Home Order","log(COVID county cases per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            notes.align="l",
            star.cutoffs = c(.1, .05,.01,.001),
            
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effects of Residing in Urban DMA in Self-Reported COVID-19 Social Distancing Behaviors",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#### Table S24

concern.la <- lm(COVIDCONCERNR1 ~ Top25Pop+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , data=datlocal)
concern.lb <- lm(COVIDCONCERNR1 ~ log(DMACoreCases)+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datlocal)
concern.lc <- lm(COVIDCONCERNR1 ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datlocal)

concern.na <- lm(COVIDCONCERNR1 ~ Top25Pop+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , data=datnolocal)
concern.nb <- lm(COVIDCONCERNR1 ~ log(DMACoreCases)+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datnolocal)
concern.nc <- lm(COVIDCONCERNR1 ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datnolocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(concern.la, concern.lb , concern.lc, concern.na, concern.nb, concern.nc,
            dep.var.caption = "",
            dep.var.labels="", 
            covariate.labels = c("Top 25 DMA","log(COVID cases in DMA County per capita)","Top 25 Covid Core", "Stay at Home Order","log(COVID cases in County)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("Local News Consumers", "Non-Local News Consumers"),
            column.separate = c(2, 2),
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="Pr(Very Worried) About Catching COVID-19",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


#### Table S25 

concern.la <- lm(COVIDCONCERN==4 ~ Top25Pop+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , data=datlocal)
concern.lb <- lm(COVIDCONCERN==4 ~ log(DMACoreCases)+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datlocal)
concern.lc <- lm(COVIDCONCERN==4 ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datlocal)

concern.na <- lm(COVIDCONCERN==4 ~ Top25Pop+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , data=datnolocal)
concern.nb <- lm(COVIDCONCERN==4 ~ log(DMACoreCases)+ stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datnolocal)
concern.nc <- lm(COVIDCONCERN==4 ~ Top25CovidCore + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED , datnolocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(concern.la, concern.lb , concern.lc, concern.na, concern.nb, concern.nc,
            dep.var.caption = "",
            dep.var.labels="", 
            covariate.labels = c("Top 25 DMA","log(COVID cases in DMA County)","Top 25 Covid Core", "Stay at Home Order","log(COVID cases in County per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("Local News Consumers", "Non-Local News Consumers"),
            column.separate = c(2, 2),
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="Pr(Very Worried) About Catching COVID-19",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)

#### Table S26


ind$localnewswatch <- NA
ind$localnewswatch[ind$MEDIA2==2] <- 1
ind$localnewswatch[ind$MEDIA2==3] <- 1
ind$localnewswatch[ind$MEDIA2==1] <- 0
ind$localnewswatch[ind$MEDIA2==4] <- 0

# PARTY7 
# 1: strong democrat, 2: not very strong dem 3: lean dem,
# 4: ind, 5: lean republican, 6: not very strong republican 7: strong republican 

concern.1 <- lm(COVIDCONCERNR1 ~ Top25Pop*localnewswatch*PARTY7 + stayathome_surveyweek + log_covidpercapitanew + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=ind)
stayhome.1 <- lm(COVIDBEHAVIOR_stayhome ~ Top25Pop*localnewswatch*PARTY7  + stayathome_surveyweek + log_covidpercapitanew  + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=ind)
maskoutside.1 <- lm(COVIDBEHAVIOR_maskoutside ~ Top25Pop*localnewswatch*PARTY7+ stayathome_surveyweek + log_covidpercapitanew  + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=ind)


table <- capture.output({ # Store the stargazer output in a string
  stargazer(concern.1, stayhome.1, maskoutside.1,  
            dep.var.caption = "",
            dep.var.labels=c("Concern", "Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("Top 25 DMA","Local News Consumer", 
                                 "Party ID (1:Strong Dem - 7:Strong Rep)", 
                                 "Stay at Home Order",
                                 "log(COVID cases in county per capita)",
                                 "Age","Weekly Church Attend","Child at Home",
                                 "Parent in Elderly Home","HS Educ. or Less",
                                 "College Educ. or More","Male","Unemployed", 
                                 "Top.25 x Local News Consumer", 
                                 "Top 25 x Party ID",
                                 "Local News Consumer x Party ID",
                                 "Top 25 x Local News Consumer x Party ID"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local", "No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            notes.align="l",
            star.cutoffs = c(.1, .05,.01,.001),
            
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="Self-Reported Concerns and Behavioral Changes",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


#### Table S27

load("tveyes_stations.RData") 

stations$fips <- stations$countyfips
stations$sinclair_final <- stations$sinclair + stations$sinclair2017

sinclairdf  <-  stations %>% select(fips, sinclair_final, dma, callsign, affiliation)

sinclairdf <- sinclairdf %>%                                        
  group_by(fips) %>%                         
  summarise_at(vars(sinclair_final),             
               list(total = sum)) 

sinclairdf$sinclair_dummy <- NA 
sinclairdf$sinclair_dummy[sinclairdf$total==0] <- 0  
sinclairdf$sinclair_dummy[sinclairdf$total > 0] <- 1  



full2 <- merge(full2, sinclairdf, by = "fips", all.x = TRUE)

## recoding so that the rest of the counties are coded as 0 (missing = 0)
full2$sinclair_new <- 0 
full2$sinclair_new[full2$sinclair_dummy==1] <- 1 

library(gmodels)
CrossTable(full2$sinclair_new, full2$Top25Pop, digits=3, max.width = 5, 
           format=c("SAS","SPSS"))


library(stargazer)

# XY Plot
plot(full2$sinclair_new,full2$SIP_Index_April1*100)

boxplot(full2$sinclair_new,full2$SIP_Index_April1*100)


travel.sinclair <- lm(SIP_Index_April1*100 ~ Top25Pop + sinclair_new + stayathome_surveyweek + 
                        log_covidpercapitanew + pop.density + med.income +
                        perc.white + perc.less.hs+ rep.perc16, data = full2)
summary(travel.sinclair)


table <- capture.output({ 
  stargazer(travel.sinclair,
            dep.var.labels=c("% Stay at Home"),
            covariate.labels = c("Top 25 DMA","Sinclair ownership","Stay at Home Order","log(COVID cases in County per capita)",
                                 "Pop. Density","Median Income","% White","% HS Educ. or less","% Trump Vote"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            no.space=TRUE,
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            notes.align="l",
            digits=1,
            title="The Effect of Sinclair Ownership on % Stay at Home",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.8\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)


#### Table S28 
load("tveyes_stations.RData") 

stations$fips <- stations$countyfips
stations$sinclair_final <- stations$sinclair + stations$sinclair2017

sinclairdf  <-  stations %>% select(fips, sinclair_final, dma, callsign, affiliation)

sinclairdf <- sinclairdf %>%                                        
  group_by(fips) %>%                         
  summarise_at(vars(sinclair_final),             
               list(total = sum)) 

sinclairdf$sinclair_dummy <- NA 
sinclairdf$sinclair_dummy[sinclairdf$total==0] <- 0  
sinclairdf$sinclair_dummy[sinclairdf$total > 0] <- 1  


full_new <- merge(ind, sinclairdf, by = "fips", all.x = TRUE)


## recoding so that the rest of the counties are coded as 0 (missing = 0)
full_new$sinclair_new <- 0 
full_new$sinclair_new[full_new$sinclair_dummy==1] <- 1 

# Subset to Local News Consumers
datlocal <- subset(full_new,full_new$MORELOCALNEWS > 1)

# Subset to National Only
datnolocal <- subset(full_new,full_new$MEDIA2 == 1 | full_new$MEDIA2 == 4)


concern.2 <- lm(COVIDCONCERNR1 ~ sinclair_new + Top25Pop +  stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
concern.4 <- lm(COVIDCONCERNR1 ~ sinclair_new + Top25Pop+  stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
stayhome.1 <- lm(COVIDBEHAVIOR_stayhome ~ sinclair_new + Top25Pop + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
stayhome.2 <- lm(COVIDBEHAVIOR_stayhome ~ sinclair_new + Top25Pop + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)
maskoutside.1 <- lm(COVIDBEHAVIOR_maskoutside ~ sinclair_new + Top25Pop + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datnolocal)
maskoutside.2 <- lm(COVIDBEHAVIOR_maskoutside ~ sinclair_new + Top25Pop + stayathome_surveyweek + log_covidpercapitanew + Dem + Rep + Age + ChurchWeekly + CHILD + OLDPARENTS + HSEduc + CollegeEduc + Male + UNEMPLOYED, data=datlocal)

table <- capture.output({ # Store the stargazer output in a string
  stargazer(concern.2, concern.4, stayhome.1,stayhome.2,maskoutside.1, maskoutside.2, 
            dep.var.caption = "",
            dep.var.labels=c("Concern", "Pr(Stay Home)","Pr(Wear Mask Outside)"), 
            covariate.labels = c("Sinclair Ownership", "Top 25 DMA","Stay at Home Order","log(COVID cases in county per capita)","Democrat","Republican","Age","Weekly Church Attend","Child at Home","Parent in Elderly Home","HS Educ. or Less","College Educ. or More","Male","Unemployed"),
            omit.stat=c("adj.rsq","LL","ser","f"),
            column.labels   = c("No Local", "Local","No Local", "Local", "No Local", "Local"),
            column.separate = c(1,1, 1, 1, 1,1),
            star.cutoffs = c(.1, .05,.01,.001),
            
            notes.align="l",
            star.char = c("+", "*", "**", "***"), 
            notes = c("+ p<0.1; * p<0.05; ** p<0.01; *** p<0.001"), 
            notes.append = F,
            no.space = TRUE,
            digits=2,
            digits.extra = 1,
            title="The Effect of Sinclair Ownership on Self-Reported Concerns and Behavioral Changes",
            align = TRUE)
})
table <- gsub("\\begin{tabular}","\\resizebox{0.9\\textwidth}{!}{\\begin{tabular}", table,fixed=T)
table <- gsub("\\end{tabular}","\\end{tabular}}", table,fixed=T)
cat(table)





